//$$ newmat.tex            Documentation file

MATRIX PACKAGE                                          30 July, 1991

Copyright (C) 1991: R B Davies and DSIR

Permission is granted to distribute but not to sell.


Introduction
------------

This is a description of an experimental matrix package. The package is
an upgrade of one I distributed to a number of people in 1990.


The package supports matrix types

Matrix                       (rectangular matrix)
UpperTriangularMatrix
LowerTriangularMatrix
DiagonalMatrix
SymmetricMatrix
RowVector                    (derived from Matrix)
ColumnVector                 (derived from Matrix).

Only one element type is supported.

The package includes the operations *, +, -, inverse, transpose,
conversion between types, submatrix, determinant, Cholesky and
Householder triangularisation.


This description begins with some general notes on matrix packages, then
a general description of my package and a list of features and problems.
This is followed by details where you can get a copy and the detailed
documentation.


What we want from a matrix package
---- -- ---- ---- - ------ -------

A matrix package needs to provide

1.   Evaluation of matrix expressions in a form familiar to
     scientists and engineers. For example  A = B * (C + D.t());
2.   Access to the elements of a matrix;
3.   Access to submatrices;
4.   Common elementary matrix functions such as determinant and trace;
5.   A platform for developing advanced matrix functions such as eigen-
     value analysis;
6.   Good efficiency and storage management.
7.   Graceful exit from errors.

It may also provide

8.   A variety of types of elements (eg real and complex);
9.   A variety of types of matrices (eg rectangular and symmetric);

It may target small matrices (say 3 x 3), or medium sized matrices, or
very large matrices.

I do not believe one can build a matrix package that will meet
everyone's requirements. In many cases if you put in one facility, you
impose overheads on everyone using the package. This both is storage
required for the program and in efficiency. Likewise a package that is
optimised towards handling large matrices is likely to become less
efficient for very small matrices where the administration time for the
matrix may become significant compared with the time to carry out the
operations.

It is better to provide a variety of packages (hopefully compatible) so
that most users can find one that meets their requirements. My package
is intended to be one of these packages; but not all of them.

I don't think it is obvious what is the best way of going about
structuring a matrix package. And I don't think you can figure this
out with "thought" experiments. Different people have to try out
different approaches. And someone else may have to figure out which is
best. Or, more likely, the ultimate packages will lift some ideas from
each of a variety of trial packages. So, I don't claim my package is an
"ultimate" package, but simply a trial of a number of ideas.

But I do hope it will meet the immediate requirements of some people who
need to carry out matrix calculations.


The present package
--- ------- -------

The package described here is intended to meet the requirements of
someone carrying out statistical analysis and related calculations.

Specifically it caters for rectangular matrices, upper and lower
triangular matrices, diagonal matrices, symmetric matrices, and row and
column vectors. Only one element type (either float or double) is
allowed. A later version may include band matrices and/or a limited
sparse matrix capability.

It is intended for matrices in the range 4 x 4 to about 80 x 80. The
upper limit is imposed by the maximum number of elements that can be
contained in a single array (8000 doubles in some machines). 

Operations provided include +, -, *, inverse, transpose, scaling by a
scalar, conversion between types, submatrix, determinant, trace, sum of
squares of elements, Householder triangularisation, Cholesky
decomposition. Singular value decomposition is to come but is not in the
present package.

The package is designed for version 2.0 of C++. It works with Sun C++,
Turbo C++, Borland C++, Glockenspiel C++ (2.00a) on a PC. There are
problems with Zortech C++ (version 2.12) and JPI C++.


Features
--------

Matrix expressions are evaluated in two stages. In the the first stage a
tree representation of the expression is formed. For example (A+B)*C is
represented by

                    *
                   / \
                  +   C
                 / \
                A   B

The expression is evaluated recursively during the second stage. The two
stage approach means that an expression can be scanned before evaluation
and the best method of evaluation found.

The package recognises combinations of the form  A.i()*B  where i() is
the inverse operation and avoids evaluating the inverse explicitly.

The package recognises and takes advantage of a matrix appearing on both
sides of an '=' as in  A=B-A; resulting in faster execution and less
temporary storage.  There is no need to provide an operator += since 
A=A+B  is automatically recognised as an "add into" operation.

Matrices generated temporarily when an expression is evaluated are
destroyed or recycled as soon as their values are no longer required.
Each matrix has a tag which indicates if it is temporary so that its
memory is available for recycling or deleting.

Further possibilities not yet included are to recognise A.t()*A and
A.t()+A as symmetric or to improve the efficiency of evaluation of
expressions like A+B+C, A*B*C, A*B.t()  [t() denotes transpose].

The package attempts to solve the problem of the large number of
versions of the binary operations required when one has a variety of
types. With n types of matrices the binary operations will each require
n-squared separate algorithms. Some reduction in the number may be
possible by carrying out conversions. However the situation rapidly
becomes impossible with more than 4 or 5 types.

In this package each matrix type includes routines for extracting
individual rows or columns. Only a single algorithm is then required for
each binary operation. The rows can be located very quickly since most
of the matrices are stored row by row. Columns must be copied and so the
access is somewhat slower. The alternative approach of using iterators
will be slower since the iterators will involve virtual functions.

In fact, I provide two algorithms for operations like + . If one is
adding two matrices of the same type then there is no need to access the
individual rows or columns and a faster general algorithm is
appropriate.

The original version of the package did not use this access by row or
column method and provided the multitude of algorithms for the
combination of different matrix types. The code file length turned out
to be just a little longer than the present one when providing the same
facilities with 5 distinct types of matrices. It would have been very
difficult to increase the number of matrix types in the original
version. Apparently 4 to 5 types is about the break even point for
switching to the approach adopted in the present package.


Problems and limitations
-------- --- -----------

The package does not have graceful exit from errors. There seems little
point in developing this until exceptions become available in C++.

The package does not have delayed copy. In most situations this is not a
significant disadvantage. It is a problem with functions returning a
matrix; see the detailed documentation.

The package does not readily handle matrices declared constant. Every
operation must have the option of recycling the memory of matrices in
its arguments. So it is not possible to declare these arguments as
constant. It is too complicated to declare versions of each operation
for constant and non-constant arguments and there doesn't seem to be a
way of doing automatic conversions. In the present package one needs to
carry out an explicit conversion.

The package is long with a large .h file. I do not recommended it if you
are short of storage and require only matrices and vectors and do not
require triangular or symmetric matrices. If you need to deal with
several types of matrices then I am pretty sure you need the row and
column operations used in this package. I am not so sure about the two
stage approach to evaluating matrix expressions.

It is also not recommended when you are dealing with very small matrices
when it becomes inefficient or very large matrices when memory
management becomes a problem (I may introduce a HugeMatrix type or alter
the memory storage method in a later version).

There is a problem with overloading a function to access a variety of
matrix types. This works fine if you access only matrices with the type
known to the compiler. But it does not work with matrix expressions
since the compiler cannot determine the type.

Matrix types of expressions are determined at run-time so the compiler
will not detect errors of the type Matrix M; DiagonalMatrix D; .... D=M;
Such errors will be detected when the statement is executed.

Symmetric matrices are not handled very efficiently (yet) since complete
rows are not stored explicitly.

It will be non-trivial to include band matrices or sparse matrices in
this package.


How to get a copy of this package
--- -- --- - ---- -- ---- -------

I am putting copies on at least some of Compuserve (Borland library, zip
format), SIMTEL20 (zip format), comp.sources.misc on Internet (shar
format), and on one of the program libraries at Victoria University,
Wellington.

-------------------------------------------------------------------------------


Detailed Documentation
======== =============

Copyright (C) 1989,1990,1991: R B Davies and DSIR

Permission is granted to use or distribute but not to sell.

Users must understand that this is a test version; there may still be
bugs and errors. Use at your own risk. Neither I nor DSIR take any
responsibility for any errors or omissions in this package or for any
misfortune that may befall you or others as a result of its use.

Please report bugs to me at

    robert@am.dsir.govt.nz

or

    Compuserve 72777,656

The matrix inverse routine is adapted from "Numerical Recipies in C" by
Press, Flannery, Teukolsky, Vetterling, published by the Cambridge
University Press.

Some other code is adapted from routines in "Handbook for Automatic
Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published
by Springer Verlag. 


The package supports matrix types

Matrix                       (rectangular matrix)
UpperTriangularMatrix
LowerTriangularMatrix
DiagonalMatrix
SymmetricMatrix
RowVector                    (derived from Matrix)
ColumnVector                 (derived from Matrix).

Only one element type is supported.

The package includes the operations *, +, -, inverse, transpose,
conversion between types, submatrix, determinant, Cholesky and
Householder triangularisation.


Customising
-----------

I use .hxx as the suffix of definition files and .cxx as the suffix of
C++ source files. This does not cause any problems with the compilers I
use except that Borland and Turbo need to be told to accept any suffix
as meaning a C++ file rather than a C file.

Each file accessing the matrix package needs to have file newmat.hxx 
#included  at the beginning. Files using matrix applications (Cholesky
decomposition, Householder triangularisation) need newmatap.hxx instead
or as well. Glockenspiel also needs to have include.hxx #included before
newmat.hxx. In addition the type of the matrix elements needs to be set
with either

   typedef double real;

or

   typedef float real;

before the include statements. If the second statement is used the
typedef statements in each of the .cxx files making up the package will
need to be updated.

In the next version of the package you also will need to edit the
#define in file precisio.hxx to reflect the use of float rather than
double. One  .cxx file will include the statement:

   typedef long double long_real;

which should be changed to

   typedef double long_real;

The include.hxx file contains details of the include files appropriate
for the compiler you are using. #define ATandT if you are using AT&T
C++. You may also need to edit this file in other ways to suit your
compiler.


Constructors
------------

To construct an m x n matrix, A, (m and n are integers) use

    Matrix A(m,n);

The UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and
DiagonalMatrix types are symmetric. To construct an n x n matrix use,
for example

    UpperTriangularMatrix U(n);

Likewise the RowVector and ColumnVector types take just one argument in
their constructors:

    RowVector RV(n);

You can also construct vectors and matrices without specifying the
dimension. For example

    Matrix A;

In this case the dimension must be set by an assignment statement or a
re-dimension statement.

You can also use a constructor to set a matrix equal to another matrix
or matrix expression.

    Matrix A = U;

    Matrix A = U * L;

Only conversions that don't lose information are supported - eg you
cannot convert an upper triangular matrix into a diagonal matrix.


Elements of matrices
-------- -- --------

Elements are accessed by expressions of the form A(i,j) where i and j
run from 1 to the appropriate dimension. Access elements of vectors with
just one argument.

This is different from the earlier version of the package in which the
subscripts ran from 0 to one less than the appropriate dimension. Use
A.element(i,j) if you want this earlier convention.


Matrix copy
------ ----

The operator = is used for copying matrices, converting matrices, or
evaluating expressions. For example

    A = B;  A = L;  A = L * U;

Only conversions that don't lose information are supported. The
dimensions of the matrix on the left hand side are adjusted to those of
the matrix or expression on the right hand side. Elements on the right
hand side which are not present on the left hand side are set to zero.

The operator << can be used in place of = where it is permissible for
information to be lost.

For example

    SymmetricMatrix S; Matrix A;
    ......
    S << A.t() * A;

is acceptable whereas

    S = A.t() * A;                            // error

will cause a runtime error since the package does not (yet) recognise
A.t()*A as symmetric.

Note that you can not use << with constructors. For example

    SymmetricMatrix S << A.t() * A;           // error

does not work.

A third copy routine is used in a similar role to =. Use

    A.Inject(D);

to copy the elements of D to the corresponding elements of A but leave
the elements of A unchanged if there is no corresponding element of D
(the = operator would set them to 0). This is useful, for example, for
setting the diagonal elements of a matrix without disturbing the rest of
the matrix.

Both << and Inject can be used with submatrix expressions on the left
hand side. See the section on submatrices.

To set the elements of a matrix to a scalar use operator =

    real r; Matrix A(m,n);
    ......
    Matrix A(m,n); A = r;

You can load the elements of a matrix from an array:

    Matrix A(3,2);
    real a[] = { 11,12,21,22,31,33 };
    A << a;

This construction cannot check that the numbers of elements match correctly.
This version of << can be used with submatrices on the left hand side.


Unary operators
----- ---------

The package supports unary operations

change sign of elements            -A
transpose                          A.t()
inverse (of square matrix A)       A.i()


Binary operations
------ ----------

The package supports binary operations

matrix addition                    A+B
matrix subtraction                 A-B
matrix multiplication              A*B
equation solve (square matrix A)   A.i()*B

In the last case the inverse is not calculated.

Notes:

If you are doing repeated multiplication. For example A*B*C, use
brackets to force the order to minimize the number of operations. If C
is a vector and A is not a vector, then it will usually reduce the
number of operations to use A*(B*C).

The package does not recognise B*A.i() as an equation solve. It is
probably better to use (A.t().i()*B.t()).t() .


Combination of a matrix and scalar
----------- -- - ------ --- ------

The following expression multiplies the elements of a matrix A by a
scalar f:  A * f; Likewise one can divide the elements of a matrix A by
a scalar f:  A / f;

The expressions  A + f and A - f add or subtract a rectangular matrix of
the same dimension as A with elements equal to f to or from the matrix
A.

In each case the matrix must be the first term in the expression.
Expressions such  f + A  or  f * A  are not recognised.


Scalar functions of matrices
------ --------- -- --------
            
    A.Nrows()
    A.Ncols()
    A.SumSquare()
    A.Trace()
    A.LogDeterminant()

have the obvious meanings. A.LogDeterminant() returns a value of type
LogandSign. If L is of type  LogAndSign  use

    L.Value()    to get the value of the determinant
    L.Sign()     to get the sign of the determinant (values 1, 0, -1)
    L.LogValue() to get the log of the absolute value.


Submatrix operations
--------- ----------

A.SubMatrix(fr,lr,fc,lc)

This selects a submatrix from A. the arguments  fr,lr,fc,lc  are the
first row, last row, first column, last column of the submatrix with the
numbering beginning at 1. This may be used in any matrix expression or
on the left hand side of << or Inject. Inject does not check no
information loss in this case. You can also use the construction

real c; .... A.SubMatrix(fr,lr,fc,lc) << c;

to set a submatrix equal to a constant.

The follwing are variants of SubMatrix:

    A.SymSubMatrix(f,l)             //   This assumes fr=fc and lr=lc.
    A.Rows(f,l)                     //   select rows
    A.Row(f)                        //   select single row
    A.Cols(f,l)                     //   select columns
    A.Col(f)                        //   select single column

In each case f and l mean the first and last row or column to be
selected.

If SubMatrix or its variant occur on the left hand side on an = or << or
within an expression its type is as follows

A.Submatrix(fr,lr,fc,lc):           Row or Column vector type of A is
                                    vector otherwise type Matrix;
A.SymSumMatrix(f,l):                Same type as A
A.Rows(f,l):                        Type Matrix
A.Row(f):                           Type RowVector
A.Cols(f,l):                        Type Matrix
A.Col(f):                           Type ColumnVector


Change dimensions
------ ----------

The following operations change the dimensions of a matrix. The values
of the elements are lost.

    A.ReDimension(nrows,ncols);     // for type Matrix
    A.ReDimension(n);               // for all other types


Change type
------ ----

The following functions interpret the elements of a matrix
(stored row by row) to be a vector or matrix of a different type. Actual
copying is usually avoided where these occur as part of a more
complicated expression.

    A.CopyToRow()
    A.CopyToColumn()
    A.CopyToDiagonal()
    A.CopyToMatrix(nrows,ncols)
    A.c()

The last expression is used in matrix expressions when A is of a const
type.


Multiple matrix solve
-------- ------ -----

If A is a square or symmetric matrix use

    CroutMatrix X = A;                // carries out LU decomposition
    Matrix AP = X.i()*P; Matrix AQ = X.i()*Q;
    LogAndSign ld = LogDeterminant(X);

rather than

    Matrix AP = A.i()*P; Matrix AQ = A.i()*Q;
    LogAndSign ld = LogDeterminant(A);

since each operation will repeat the LU decompostion.


Memory management
------ ----------

The package does not support delayed copy. Several strategies are
required to prevent unnecessary matrix copies.

Where a matrix is called as a function argument use a constant
reference. For example

    YourFunction(const Matrix& A)

rather than

   YourFunction(Matrix A)

Constant matrices cannot be used in matrix expressions so if you wish to
use A in an expression within this function use A.c() rather than A.

A second place where it is desirable to avoid unnecessary copies is when
a function is returning a matrix. When returning a matrix preceed the
return statement with a Release statement (or a ReleaseAndDelete
statement if the matrix was created with new). For example

   Matrix A;
   ......
   A.Release(); return A;

or

   Matrix* m = new Matrix;
   ......
   m->ReleaseAndDelete(); return *m;

Note that some compilers (Zortech, for example) do not use the
initialise operation when returning an object. This construction must
not be used in these cases.

The memory management methods used in this package often do a better job
than delayed copy for minimising unnecessary storage of matrices. I saw
delayed copy as an additional complexity that was probably not worth
including.


Cholesky decomposition
-------- -------------

Suppose S is Symmetric and positive definite. Then there exists a unique
lower triangular matrix L such that L * L.t() = S. To calculate this use

SymmetricMatrix S;
......
LowerTriangularMatrix L = Cholesky(S);

You will need to delete the Release() statement from the code of
Cholesky if you are using a compiler that avoids the initialise
instruction when returning an object.


Householder triangularisation
----------- -----------------

Start with matrix

       / X    0 \      s
       \ Y    0 /      t

         n    s

The Householder triangularisation post multiplies by an orthogonal
matrix Q such that the matrix becomes

       / 0    L \      s
       \ Z    M /      t

         n    s

where L is lower triangular. Note that X is the transpose of the matrix
sometimes considered in this context.

This is good for solving least squares problems: choose b (matrix or row
vector) to minimize the sum of the squares of the elements of

         Y - b*X

Then choose b = M * L.i();

Two routines are provided:

    HHDecompose(X, L);

replaces X by orthogonal columns and forms L.

    HHDecompose(X, Y, M);

uses X from the first routine, replaces Y by Z and forms M.


